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Abstract. Within the recently proposed doped-carrier representation of the 
projected lattice electron operators we derive a full Ising version of the i-J model. 
This model possesses the global discrete Z2 symmetry as a maximal spin symmetry of 
the Hamiltonian at any values of the coupling constants, t and J. In contrast, in the 
spin anisotropic limit of the t-J model, usually referred to as the t-Jz model, the global 
SU{2) invariance is fully restored at = 0, so that only the spin-spin interaction has 
in that model the true Ising form. We discuss a relationship between those two models 
and the standard isotropic t-J model. We show that the low-energy quasiparticles 
in all three models share the qualitatively similar properties at low doping and small 
values of J/ 1. The main advantage of the proposed Ising t-J model over the t-Jz one 
is that the former allows for the unbiased Monte Carlo calculations on large clusters 
of up to 10'^ sites. Within this model we discuss in detail the destruction of the 
antiferromagnetic (AF) order by doping as well as the interplay between the AF order 
and hole mobility. We also discuss the effect of the exchange interaction and that of the 
next nearest neighbour hoppings on the destruction of the AF order at finite doping. 
We show that the short-range AF order is observed in a wide range of temperatures 
and dopings, much beyond the boundaries of the AF phase. We explicitly demonstrate 
that the local no double occupancy constraint plays the dominant role in destroying 
the magnetic order at finite doping. Finally, a role of inhomogeneities is discussed. 
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1. Introduction 

Since tlie discovery of the high-temperature superconductors (HTSC), many theoretical 
investigations have been focused on the study of the doping evolution from the 
antiferromagnetically ordered Mott insulator to the BCS-type superconductor. It is 
commonly assumed that this evolution is the key ingredient for the understanding of the 
physics behind high-temperature superconductivity and it can adequately be addressed 
in terms of the two-dimensional {2d) t-J lattice model. This model is believed to capture 
the essential low-energy physics of doped Mott insulators driven by strong electron 
correlations. Although, it is generally accepted that strong electron correlations play 
an important role close to half filling, it remains unclear to what extent HTSC can be 
described entirely in terms of this minimal electronic model. 

One of the reasons for this lack of clarity is that, despite its simplicity, the exact 
properties of the t — J model, apart from a few limiting cases, are still unknown. The 
vast majority of the results away from half filling have been obtained for one or two holes 
introduced into the AF background. This problem has been thoroughly analyzed with 
the help of various analytical and numerical approaches [Ill2l[3],[5],|4l[6l[7l[8],[9l [TOl [TT]. 
Results for larger dopings are less comprehensive. Here, one of the important problems 
concerns the robustness of the AF order against doping. Most of the theoretical 
approaches predict, that the long range AF order persists up to much larger dopings 
than observed in cuprates [12l [13] . It has recently been suggested that including into 
consideration the nearest neighbour hopping may help in resolving this discrepancy |14j . 
However, a reliable and well controlled analytical treatment of the isotropic t-J model 
poses a severe technical problem: it is very hard to analytically deal with the local 
no double occupancy (NDO) constraint. On the other hand, because of this constraint, 
numerical treatment is available only on rather small lattice clusters, which immediately 
rises the problem of the finite-size effects and, consequently, of the thermodynamical 
stabihty of the results obtained. 

An interesting question then arises: is there available a modification of the isotropic 
t — J model that allows for unbiased numerical treatment and at the same time captures 
the essential properties of the original t — J model at least for a certain parameter range? 
Some of the computational difficulties related with the t-J model may be overcome by 
means of investigations of its anisotropic limit, i.e., the t-Jz model. Numerical work 
has suggested that the t — J and t — J^ models have many similar properties [15]. In 
particular, the stripes, pairing, and Nagaoka states found in the t — J^ model are very 
similar to those of the isotropic t — J model [16], [171 [HI [19], [20]. In general, despite some 
significant differences between both the models [Sni [21] , at the low-energy scale of order 
J (-C t) it is reasonable to consider the background spin configuration to be frozen with 
respect to the hole dynamics time scale. In this case the properties of the low-energy 
quasiparticle excitations in the t — J model are at least qualitatively similar to those 
in the anisotropic t — Jz model. Although this model is more amenable to numerical 
calculations, again only rather small lattice clusters are allowed. 
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One may hope that the full Ising version of the t — J model in which the t-term 
possesses the global discrete Z2 spin symmetry rather than the SU(2) one results in a 
more tractable though still nontrivial model. However, it is not clear how such a model 
can be derived directly in terms of the Gutzwiller projected lattice electron operators, 
since those operators transform themselves in a fundamental representation of SU(2). In 
this paper we use the recently proposed doped-particle representation of the projected 
electron operators to derive the full Ising version of the t — J model. We believe that the 
results of our study capture some essential low-energy properties of the isotropic t — J 
model since the strongly correlated nature of the problem is preserved. 

The doped-particle representation of the t-J model is especially suited for 
investigations of the underdoped regime [22|, [231 [2l] . Some crucial points of this approach 
are recalled at the beginning of the following section. Although this is a slave-particle 
formulation, it differs from other similar approaches in that the NDO constraint is 
identically fulfilled at half-filling. Since the holes are the only charge degrees of freedom 
present in the system, the number of charge carriers is very small close to half-filling. 
Therefore, this approach is particularly useful for the description of the underdoped 
regime, whereas its applicability to the description of strongly overdoped cuprates is 
much more involved. 

Until now, the doped particle representation has been analyzed only within mean- 
field approximations [221 EH]. Our aim is to go beyond this approximation in such a 
way, that the slave-particle constraint is exactly taken into account. We demonstrate 
that this can be achieved for very large clusters of the order of 10^ lattice sites, provided 
the SU(2) symmetry is broken down to Z2 not only in the spin-spin interaction term 
but also in the hopping term. Since the effectiveness of the resulting calculations is 
independent of the translational invariance, this approach also allows one to investigate 
the role of inhomogeneities which are expected to play an important role in the cuprate 
compounds [25] . 

The paper is organized as follows. Section 2 describes in detail the derivation of the 
theoretical model. Section 3 comprises the results of the numerical calculations. Section 
4 lists our conclusions. 

2. Model and Approach 

2.1. Doped- carrier formulation of the t-J model 

We start with the t-J Hamiltonian on a square lattice [26] 

Htj = -Y. tiMaCja + J Y^^QiQj - (1) 

ija- (ij) 

where Cia = PciaP = Cto-(l — ni-a) is the projected electron operator (to exclude the 
on-site double occupancy), Qi = J2a,a' cl^'''aa'Cia' , = 3/4, is the electron spin operator 
and fii = PuiP = Ui-^ + rii^ — 2nj|njj^. Hamiltonian ([1]) contains a kinetic term with the 
hopping integrals tij and a potential J describing the strength of the nearest neighbour 
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spin exchange interaction. At every lattice site the Gutzwiller projection operator 
P = Hii^ ^ nia-ni_^) projects out the doubly occupied states | ti) thereby reducing 
the quantum Hilbert space to a product of 3- dimensional spaces spanned by the states 
|0)j, I |)j and I Physically this modification of the original Hilbert space results in 
strong electron correlation effects. The crucial local no double occupancy constraint 
is rigorously incorporated into Eq. ([1]). However, this is achieved at the expense of 
introducing the constrained electron operators, c|^, that obey much more complicated 
commutation relations than the conventional "unconstrained" fermion operators. It 
should be stressed that it is precisely close to half filling where the Gutzwiller projection 
is of a crucial importance: the projected electron operator c}^ in this regime significantly 
differs from the bare electron operator c}^ (right at half filling cj^ = 0). 

A natural question then arises: is it possible to rewrite the t-J Hamiltonian in terms 
of the conventional fermion and spin operators in such a way that the NDO constraint 
for the lattice electrons transforms itself into the one that can be treated, close to half 
filling, in a controlled way? Recently, it has been shown that the t-J Hamiltonian can 
indeed be represented in that form [23l [22| [23]. 

For the reader's convenience we sketch below the main points of this scheme. The 
basic idea behind this approach is to assign fermion operators to doped carriers (holes, 
for example) rather than to the lattice electrons. The t-J Hamiltonian is expressed then 
in terms of the lattice spin operators. Si, and doped carrier operators represented by 
spin-1/2 charged fermions, di„. 

To accommodate these new operators one obviously needs to enlarge the original 
onsite Hilbert space of quantum states. This enlarged space is characterized by the state 
vectors \aa) with a =il, -l| labeling the spin projection of the lattice spins and a = 0, j, | 
labeling the dopon states (double occupancy is not allowed). In this way the enlarged 
Hilbert space becomes 

nr' = {\t o)„ 1 4 o)„ I I ^T)., I ^T)., I ^i).}, (2) 

while in the original Hilbert space we can either have one electron with spin a or a 
vacancy: 

^ = {|T).,li).,|0),}. (3) 
The following mapping between the two spaces is then defined [22] : 

IT). --|^0)„ |i).^UO)„ (4) 

|0), - (5) 

The remaining states in the enlarged Hilbert space, (| + | JJ-t)*) /v^? I ltT)j! I 
are unphysical and should therefore be removed in actual calculations. In this mapping, 
a vacancy in the electronic system corresponds to a singlet pair of a lattice spin and a 
dopon, whereas the presence of an electron is related to the absence of a dopon. 
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The spin-dopon representation of the t-J Hamiltonian then reads 

ija (ij) 

with di(j = di^{l — d\_^di^^ij) being a projected dopon operator and hf = Y, dladia- 

The apphcation of in this form should be accompanied by the implementation 
of the constraint to eliminate the unphysical states, 

5,M, + -^fit = 0, (7) 

where = J2a,a' dl^Taa' dia' stands for the dopon spin operator so that Qi = Si + Mi. 
Note the important factor of 2 in front of the first term in Eq. ([6]). It originates from 
the fact that the vacancies are represented in this theory by the spin-dopon singlets 
given by Eq. ([5]). The projected lattice electron operators can be explicitly expressed 
in terms of the projected dopon operators. For example, 



^Yy2^ Sl'j dil- S^di-i , (8) 

where Pf'' = l-{SiMi + \nf) is the projection operator which eliminates the unphysical 
states from the i^^ site. 

The above representation of the t-J Hamiltonian is particularly useful for the 
description of strongly underdoped cuprates. Close to half filling nf <^ 1, so that one 
can safely drop the tilde sign off the projected dopon operators. This is due to the fact 
that in the low doping regime the probability for the realization of a doubly occupied 
dopon state is indeed very low. Despite that, the NDO constraint ([7]) must be imposed 
to eliminate the unphysical degrees of freedom that are present in this formalism at 
any finite doping. Note, however, that at half-filling the left hand side of Eq. ([7]) 
vanishes, and thus, in contrast to the original NDO constraint for the lattice electrons, 
this equation turns into a trivial identity||| Additionally, in this regime one can neglect 
both the hole-hole interactions represented by the MiMj term, and the hfh'^ couplings. 
Note also that the insulating phase in this representation is directly associated with the 
absence of charged particles. 



2.2. The t-Jz and the Ising t — J models 

So far, the doped-particle representation of the t-J model has been analyzed only within 
the mean-field approximations [23l [22], [2l]. Our aim is to go beyond the mean-field 
analysis and to carry out calculations for large enough systems, to make sure that the 
finite-size effects are truly negligible. Let us start with the anisotropic limiting case of 
the t-J model with the spins polarized only along the z-component, namely, the t-J^ 

X The original local NDO constraint for the lattice electrons, cL'^ic — 1' I'ight at half filling reads 
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model. This model is of interest in itself since it captures some essential physics of 
strong electron correlations. 

The t-Jz model can be considered as a limiting case of the t-J model ([I]) which 
has an Ising rather than a Heisenberg spin interaction: 

Ht-j-. = - E ^ijcLcja + JzY. (QiQj - ' (9) 

ijcr (ij) 

Here represent the Gutzwiller-projected electron operators. The global continuous 
spin SU(2) symmetry of the t-J model now reduces to the global discrete Z2 symmetry 
of the t-Jz model. Although Q^Qj interaction possesses discrete Z2 symmetry, the 
original SU(2) symmetry of all other terms of the Hamiltonian is preserved. Therefore, 
the symmetry of the t-Jz model depends on whether Jz is zero or finite. Namely, for 
Jz = the SU(2) symmetry is restored again. In contrast, in the full t — J model both 
the t- and J-terms possess the same SU(2) symmetry. 

It is therefore natural to seek for a representation of the full Ising version of the 
t — J model in which the symmetry of the model does not depend on the values of 
the model parameters. Such a representation can straightforwardly be derived within 
the dopant-particle formulation of the t-J model. The physical consequences as well 
as computational advantages of such an approach will be discussed in the subsequent 
sections. 

To proceed, we start right from the original t-J Hamiltonian described in terms of 
the lattice electrons given by Eq. ([T]), where we now put Qf = = at the operator 
level. We then have for the physical electron projected operators: 

= Vtdl^f = il + S^dl 
c.i = Vfdl^Vf = {\- S!)dl (10) 

where the projection operator now reads vf=l- (25/ Mf + nf/2). It can easily be 
checked that 

Qt = (Q^y = clrQi = 0. 

The kinetic t-term built out of the physical electron operators given by Eqs. ffTOl) 
possesses the global Z2 symmetry rather than the SU(2) one. 

Accordingly, the underlying onsite Hilbert space rearranges itself in the following 
way. The operators Cj|, cj^ act on the Hilbert space TCi = {| Jj-, 0), | Jj-)T)}- These 
operators do not mix up any other states. Operator Cj^ destroys the spin-down electron 
and creates a vacancy. This vacancy is described by the state | JJ-jT)- The similar 
consideration holds for the Q| operators. Now, however, the vacancy is described by 
the state | lt, i). Those two vacancy states are related by the Z2 transformation. The 
operator (Q|)^ = |(1 — nf) produces zero upon acting on the both. The physical 
Hilbert state is therefore a direct sum Tiph = 7i| © 7i|. Under the Z2 transformation 
(t^i, S- —S-) we get ^ Hi, which results in Hph — > Hph- 

§ In the isotropic t-J model these two 2d spaces merge into a 3d SU(2) invariant physical space, where 
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As a result, one arrives at the representation ^ in which, however, the electron 
projection operators are given by Eqs. ffTOj) . All the parts of this Hamiltonian possess 
the global discrete Z2 symmetry whereas the global continuous SU(2) symmetry is 
completely lost. Close to half-filling this Hamiltonian reduces to the form. 



hI!!j^ = Ht-,j\z^ - ^tijd\„dj„ + 

ija (ij) 

+s^m; + s;mi], (ii) 

which should be accompanied by the constraint, 

2S!Ml + \ni = 0. (12) 

In order to emphasize the difference between the t-J^ and the present model we dub 
the latter the Ising t — J or for short the t — J\z2 model. 

The factor of 2 which is presented in the hopping term of the isotropic t-J 
Hamiltonian ([6]) drops out from the hopping term in the Ising t — J Hamiltonian. This 
occurs because of the fact that the equations 

4 = V2V!%,Vf\ vf = i- + \hi 

which are valid for both the t-J and t-J^ models, are in the t — J\z2 model replaced by 
the following ones. 



In practical calculations, the NDO constraint (fT2|) can be taken into account with 
the help of a Lagrange multiplier. In order to do that we introduce an additional term 
to the Hamiltonian, 



1 



\ + St) d\^d,^+{^--St]4^d,y 



(13) 



Notice that the operator (| + S^dl-^di-^ + (| ~ ^D'^li'^ii produces eigenvalues and 1, 
when acting on the onsite physical and unphysical states, respectively. Because of this, 
the global Lagrange multiplier A — > 00 enforces the NDO constraint locally. The double 
dopon occupancy of an arbitrary site results in an appearance of an unphysical state 
and hence enhances the total energy by A. Therefore, in the large-A limit all unphysical 
states are automatically eliminated and we can in this limit safely remove the tilde sign 
off the d operators. In the following section we show that this constraint is of crucial 
importance for the description of the AF order at finite doping. 

the vacancy is just an antisymmetric linear combination given by the SU(2) spin singlet (5). The 
symmetric combination splits off, since it represents an unphysical spin-triplet state. 
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2.3. Monte Carlo calculations 



The total Hamiltonian takes the form 



Hf + Hi 



{ij) 



const, 



(14) 



with 



i 



X 



Hi 



^3 



A 



QZ 



J 



(15) 



J 



(16) 



where {j)i denotes neighbouring sites of a given site i. We have neglected the hole-hole 
interaction in H, 



t-J\zo 



which is perfectly justified in the low doping regime. In order to 
verify this approximation we have carried out additional calculations with the hole-hole 
interaction being taken into account in the mean-field approximation. The difference 
is negligible and therefore we do not present them here. Note, that the interaction 
strength in H^_j^^ is exactly the same as in the standard formulation of the t-J model. 
Absence of any renormalization of the model parameters originates from the fact that 
the projection procedure is explicitly built in Eq. ( |T4l) . provided A — oo. Despite its 
complexity, with the Monte Carlo (MC) method, one can investigate the Hamiltonian 
( HM for very large systems without any approximation. Since [5*/, H^_j^^ ] = the spin 
degrees of freedom can be analyzed within the classical Metropolis algorithm. However, 
since the effective Hamiltonian f|T^ includes both fermionic as well as classical degrees 
of freedom, this algorithm needs to be modified. The procedure is as follows: 

(i) an initial configuration of {S^} is generated; 

(ii) the Hamiltonians (ITSl [T6l) are diagonalized and the free energy JF of the fermionic 
subsystem in the canonical ensemble is determined; 

two sites with opposite spins are randomly chosen; then, both the spins are 
flipped; 

step (ii) is repeated, determining new value of the free energy JF'; 

(v) if JF' < JF or exp [(JF — JF') / kT] > x, where x is a random number from the interval 
[0; 1), the new {S^} configuration is accepted, added to the ensemble and the 
procedure goes to step (iii), otherwise it goes directly to step (iii). 

It is the Metropolis algorithm, but with the internal energy in statistical weights 
replaced by the free energy of the fermionic subsystem. A detailed description of this 
approach can be found in Ref. [27]. Concurrently with the MC simulation an iterate 
procedure calculating the distribution of holes is carried out in a self-consistent way. 



m 



IV 
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Most of the numerical results have been obtained for a 20 x 20 systems with periodic 
boundary conditions. However, in order to check the influence of finite-size effects we 
have carried out calculations on clusters of up to 1600 lattice sites and with averaging 
over the boundary conditions [28] . This problem is discussed at the end of the next 
section. 

In order to eliminate the unphysical states in the Monte Carlo simulations, we have 
taken A = lOOt. Therefore, the Lagrange multiplier is by far the largest energy scale in 
the system, which guaranties the single occupancy of each lattice site. The simulations 
have been carried out in the canonical ensemble, what allows for accurate control of 
the doping level. We assume an absence of the ferromagnetic order. Namely, we take 
I]j Sf = J2i = 0, what is reflected in the third point of the MC procedure. 

The aim of the simulations is to determine how the antiferromagnetic order and 
the spectral properties are affected by doping. We start our discussion with the doping 
dependence of the spin-spin correlation function for the projected physical electron 
operators 

9ir) = ^.EEiiS! + MDis; + m;)) 

i j 

X exp [iK ■ {Ri - Rj)] 6{r -\Ri-Rj\), (17) 
where K = (tt, vr) and 

r/ X f 1 if Ixl < 0.5a, 
o(x] = < , 

[ otherwise, 

with a being the lattice constant. (...) in Eq. (fT71) means an average over the 
spin configurations generated in MC run. g{r) allows one to distinguish between 
the long range order (LRO), when it remains finite for arbitrary r, quasi long range 
order (QLRO), when g{r) decays algebraically, and a short range order (SRO), when 
g{r) decays exponentially. Calculations of the spin-spin correlation function will be 
accompanied by results obtained for a static spin-structure factor, defined as 

siQ) = ]^ E e^'^(^-^^)((5,^ + M^is; + m;)). (is) 

The third quantity that we use in the following discussion, is the hole spectral function 
given by 

A{k,uj) = Im G (fc, + zO+) , (19) 



TT 



where 



C (fc, ^) = E E exp {^k [R, - R,)} 

-1 



X {g^{Ri,Rj,z) 



), (20) 



with s(t) = 1 and s{[) = —1. Here, similarly to Eq. flTTl) . (...) indicates averaging over 
spin configurations and 

g^{R,,R,,z) = {[z-H,]-'}^, (21) 
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is the real-space Green function for a given spin configuration {Sf}. The presence 
of factors \ — s{a)Sl in Eq. (120|) follows from Eq. ffTOl) . Note that the spin-spin 
correlation function, the spin-structure factor and the spectral function are defined for 
physical electron operators, Cj. 

3. Numerical results 

3.1. Homogeneous systems 

As discussed in the preceding sections the derived representation differs from the 
standard t-J^ Hamiltonian in that the SU(2) symmetry is broken also for J = 0. In 
order to visualize the physical consequences of this difference we start with calculations 
for the one-hole case. In this regime large clusters have been analyzed numerically both 
for the t-Jz and t-J models. 




Figure 1. One-hole energy as a function of J/t for a 8 x 8 cluster at T = 0. The 
data labeled ast—J^ and t — J have been taken from Refs. [29] and [11] , respectively. 
In the latter case, it is the energy of a hole with momentum (7r/2,7r/2). The inset 
shows results obtained for 4 x 4 t— J\z2 and t-J^ systems. Here, the dashed line shows 
results presented in Ref. [30j for the t~Jz model. Since for J ^ the ferromagnetic 
order sets in, the constraint = is now relaxed. 

In the main panel of Fig. [1] we show the one-hole energy calculated at T = for 
a 8 X 8 cluster without the ^nifij term. Here, we compare our data with exact results 
obtained for a 50-site t-J^ cluster [21] as well as with recent exact results for a bulk t-J 
system [11]. In the inset of Fig. [T]we compare exact results obtained for 4 x 4 t — J\z2 
and t-Jz clusters [3D] for a wider range of the J/t ratio. In two limiting cases J/t ^ 
and J/t ^ oo the one-hole energy obtained in our approach is the same as in the t-Jz 
model. It can be explained in the following way. For J = 0, the ferromagnetic Nagaoka 
state becomes a ground state in both the approaches. In this case the propagation of a 
hole is not perturbed by the magnetic order and the one-hole energy equals —At. The 
main difference between t-Jz and t — J\z2 approaches consists in the symmetry of the 
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hopping term. Therefore, they merge in the case J/t ^ 1, when the system properties 
are determined predominantly by the same spin-spin interaction. In the regime of 
intermediate J the differences are most pronounced. The one-hole energy obtained for 
the t-J model in this regime is in between the results obtained for t-J^ and t — J\z2 
approaches. Here, the differences between t-J and t-Jz models are comparable to those 
between t-J and t — J\z2 ones. 
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Figure 2. Panels A) and B) show g{r) calculated for kT = O.lt. J = 0.2t (A) 
and J = OAt (B). The curves from the top to the bottom have been obtained for 
S = 0.02, 0.04, 0.06, 0.08. Sohd (dashed) hnes show results obtained for t' = <" = 
(t' = —0.27t and t" = 0.2t). Panel C) shows doping dependence of the static spin- 
structure factor ^(Tr, tt) at kT = OAt for t' = -0.27< and t" = 0.2t. 



Next, we investigate how the doping affects the antiferromagnetic order. The 
previous studies of the t-J model clearly indicate that the long range hopping amplitudes 
significantly modify the bandwidth and the dispersion of the quasiparticles |31j . 
Recent Green's function Monte Carlo calculations demonstrate that the next nearest 
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neighbour hopping reduces the critical doping at which the AF LRO disappears |14j . 
The importance of these results follows from the fact that in the t-J model with 
only the nearest neighbour hopping, antiferromagnetic correlations persist up to hole 
concentrations much larger than the ones observed in HTSC materials. One may expect 
that in the absence of the transverse spin-spin interaction the robustness of LRO should 
be even more pronounced. Moreover, the previous analysis of the one- and two-hole 
spectra in the t-J^ model [20] has shown that for t/ J < 5, half of the one-hole band 
width does not exceed 0.08t. Therefore, the intra-sublattice hopping should be a source 
of an important contribution to the kinetic energy even for small values of t' and t" . 
The significance of the long range hopping for the AF order is demonstrated in Fig. 
[2l where we compare the spin-spin correlation functions calculated with and without 
t' and t" . In all the figures showing g{r) we use logarithmic scale for the vertical axis. 
Therefore, for LRO, QLRO and SRO, g{r) should be represented asymptotically by a 
constant function, logarithmic function and a straight line, respectively. In the following 
5 denotes the average concentration of holes. 

One can see that the influence of the long range hopping depends, even qualitatively, 
on the value of the exchange coupling J. For a small value of J, the AF order is enhanced 
when hoppings to second and third nearest neighbours are allowed (for J = 0.2t see panel 
A) in Fig. [2]). On the other hand, for bigger values of J, these hoppings reduce the AF 
order (for J = OAt see panel B) in Fig. [2]). Such a behaviour could be explained as 
follows. In the presence of only nearest neighbour hopping there is a strong competition 
between the energy of spin-spin interaction and the hole kinetic energy. It results from 
the fact, that in this case only inter sublattice hopping is allowed. From Eq. (fT3ll 
one can then infer that it is possible only in regions where the AF order is absent. 
Then, nonzero t' and t" allow for intra sublattice hopping, thereby leading to gaining 
of the kinetic energy without destroying the AF order. This mechanism is effective for 
J < 0.25t. On the other hand, for t' = t" = and large J, holes are almost localized 
and, therefore, only weakly frustrate the AF state. The intra-sublattice hopping allows 
for the propagation of holes, which effectively reduces the AF LRO. 

The doping induced destruction of the AF LRO can directly be seen in Fig. 
|2]G, where we present spin-structure factor obtained for t' = — 0.27t and t" = 0.2t. 
This quantity is important in that it is directly accessible in, e.g., neutron scattering 
experiments. The maximal doping for which the AF state still exists strongly depends 
on the magnitude of the exchange interaction. This result contrasts with the recently 
reported Green's function Monte Carlo study of the t — J model [H], where the AF 
LRO vanishes at 5 = 0.1 and 6 = 0.13 for J = 0.2t and J = OAt, respectively. In 
our approach the experimental data for the critical doping in HTSC can be reproduced 
provided J < 0.2t. 

In order to illustrate the interplay between the AF order and the mobility of holes 
we have calculated the hole spectral functions A{k,u) (see Figs. |3] and Hj). For 
t' = t" = and small 6 one can see almost localized particles with very small dispersion. 
Similar situation occurs in the t-J^ model, but it is not the case for the t — J one, where 
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Figure 3. Spectral functions A{k,uj) calculated for kT = O.lt along the main 
symmetry lines of the Brillouin zone with <' = t" = 0. 

the the spin-flip term can undo the defects generated by a moving hole and hence it 
allows for its much higher mobility. In the present approach holes become mobile when 
doping increases, i.e., when the AF background disappears. The onset of mobile holes 
is accomplished through a gradual transfer of the spectral weight from the vicinity of 
the almost localized level. Close to half-filling the most significant transfer takes place 
in states with k = (0,0) and k = (tt, vr). However, even for relatively large doping the 
spectral functions remain broad for all the momenta. 

For non-zero t' and t" there are mobile holes even for small doping, but the spectral 
functions still remain very broad. Also in this case, doping is responsible for significant 
modification of the dispersion relation of holes. In Fig. |l]we compare A{k, u) calculated 
for 5 = 0.02, ..,0.14 with t' = -0.27t and t" = 0.2t. Along with the destruction of 
the AF LRO, there is an increasing contribution of nearest neighbour hopping to the 
hole kinetic energy. For 6 = 0.02 the peaks in spectral functions can be fitted by 
the dispersion relation with t = 0, whereas for 6 = 0.14 the AF correlations hardly 
influence the nearest neighbour hopping. Note that such a substantial modification of 
the dispersion relation may change the topology of the Fermi surface. Doping affects 
not only the effective dispersion relation, but also frustrates the AF background. The 
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latter effect is responsible for strong broadening of the spectral functions that is visible 
in Fig. m Comparison of Fig. [2] with Figs. [3] and H] demonstrates that the mobility of 
holes and destruction of the AF LRO are mutually connected with each other in the 
sense that mobility affects AF LRO and, vice versa, AF order affects the mobility of 
holes. 

We now turn our attention to the temperature dependence of the spin-spin 
correlation function and the spectral properties of holes. It is known that the 
Neel temperature drops rapidly when the parent compounds of the high-temperature 
superconductors are doped with holes. Similar behaviour can be inferred from Fig. [5l 
where temperature dependence of g{r) is presented for different doping levels. One can 
note that doping strongly reduces the LRO, whereas its influence on the SRO is much 
weaker. In particular, the nearest neighbour correlation functions g{l) calculated for 
5 = 0.02 and 6 = 0.06 are qualitatively and quantitatively close to each other. These 
results suggest that, the AF SRO should be observed in a wide range of temperatures 
and dopings, much beyond the boundaries of the AF phase. It remains in agreement with 
recent experiments on high-temperature superconductors suggesting that with doping, 
the long-range Neel order gives way to short-range order with a progressively shorter 
correlation length. As a result, at optimal doping the static spin correlation length is 
no more than two or three lattice spacing |32j . 

The discussed above correlation function and the spin-structure factor describe 
the background composed of the localized spins, which, as mentioned in the preceding 
paragraphs, is up to some degree affected by the motion of doped holes. Therefore, the 
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Figure 5. g{r) as a function of temperature for J = Q.2t and 5 = 0.02 (upper panel) 
and S = 0.06 (lower panel). The lines from the top to the bottom show g{l), g{2) and 
.g(10). t' = -G.27t and t" = 0.2t have been assumed. 



reduction of the antiferromagnetic correlations has to be observed also in the dynamics 
of the carriers. It is shown in Fig. [6l where we demonstrate the temperature dependence 
of the spectral functions. When the temperature increases, the number of spin defects 
in the Neel state increases as well, and this enables the nearest neighbour hopping, 
thereby allowing holes to lower their kinetic energy. This mechanism leads to trapping 
of holes in the regions of broken antiferromagnetic bonds and forming ferromagnetic 
spin polarons, where the hole hopping does not frustrate the spin background. The 
contribution of the nearest neighbour hopping becomes visible in the spectral functions, 
where the increase of the temperature causes a significant broadening of the spectral 
lines. Similarly to the spectral functions obtained for a single hole in the t-t'-t"-J model 
[33], the width of the peaks in the A{k,uj) is too small, when compared to the results 
of the angle-resolved spectroscopy (ARPES) measurements [3lj on Sr2Cu02Cl2. It has 
recently been argued, that strong electron-phonon interaction [351 EH] may explain the 
very broad peaks observed in the insulating copper oxides [371 EH] • 

In the present approach, holes interact with spins through the intersite interaction 
of strength J as well as through the onsite constraint with the Lagrange multiplier A, 
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and both of these interactions can be responsible for the destruction of the AF LRO. 
In order to determine the underlying mechanism we have carried out simulations taking 
into account one of these interactions at a time. The results shown in the upper panel 
of Fig. [7] clearly demonstrate that the constraint plays the dominating role in the 
destruction of the AF state. When both these interactions are taken into account the 
LRO is completely destroyed and only AF SRO can be observed for 5 = 0.08. However, 
if we ignore the constraint LRO is restored. The results obtained without the intersite 
spin-hole interaction (S'f MJ terms) are almost indistinguishable from the ones obtained 
with both the interactions. Note, that for A = the NDO constraint is completely 
neglected. This illustrates the importance of the constraint for a realistic description 
of the AF order at finite doping. Additionally, in order to determine the value of the 
Lagrange multiplier for which the results converge and the NDO constraint is fulfilled 
we have calculated g{r) for a wide range of A. The results are presented in the lower 
panel of Fig. [3 One can see that the assumed value A = lOOt is large enough to enforce 
the constraint. After explaining the role of the farther neighbour hopping we restrict 
the following analysis to the case of f = — 0.27t and t" = 0.2t. 
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Figure 7. g{r) for J = 0.2t, /cT = O.lt, (5 = 0.08, t' = -0.27t and t" = 0.2t. In the 
upper panel the topmost curve has been calculated with A = 0. The other curves show 
results obtained for A ~ lOOt with and without the spin-hole exchange interaction. The 
lower panel shows g{r) for different values of A. Note, that the results for A = lOOt 
and A = 300t are almost indistinguishable. 



3.2. Inhomogeneous systems 

Since our method works for systems with broken translational invariance, it is tempting 
to apply it as well to inhomogeneous systems. It has recently been shown with the help of 
scanning tunneling spectroscopy that nanoscale electronic inhomogeneity is an inherent 
feature of many groups of high-temperature superconductors. By a direct probing of the 
local density of states, these methods reveal strong spatial modulation of the energy gap 
in the superconducting Bi-based compounds [39] . Very recently, STM experiments have 
shown a strong correlation between position of the dopant atoms and all manifestations 
of the nanoscale electronic disorder [iQlHT]. Thus, these experiments proved essentially 
that the impurities were the source of the inhomogeneity. On the other hand, they 
revealed a very important feature: there is a positive correlation between the magnitude 
of a gap and the position of an out-of-plane oxygen atoms [IQ]. These are the atoms 
which have been doped into the insulating parent compound in order to introduce holes 
to the Cu02 planes. The gap-impurity correlation has been explained as a result of the 
inhomogeneity-enhanced exchange interaction in the t-J model [42]. Assuming that 
purely electronic models contain the essential physics of cuprates, the same interaction 
is responsible for both superconductivity and the AF order. Therefore, inhomogeneity 
may affect the AF state as well. Additionally, localization of holes by the electrostatic 
potential of out-of-plane oxygen atoms may also affect the AF order, since the hopping 
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of holes frustrates the AF LRO. In the following, we investigate the role played by these 
mechanisms in the strongly underdoped regime. In order to investigate the latter one, we 
extend the Hamiltonian (|T4l) by adding a term responsible for inhomogeneity-induced 
diagonal disorder 

Hl.j\,, Hlj^,^ + ^^dld^a, (22) 

ia 

where 

{V, if there is an out-of-plane oxygen atom 
above site i, 
0, otherwise. 

Since in HTSC each doped oxygen atom introduces one hole in the Cu02 plane, we 
have carried out calculations for the number of impurities equal to the number of holes. 
Technically, for each Monte Carlo simulation we generate a random configuration of the 
out-of-plane oxygen atoms and keep it frozen during the whole run. In that way both 
holes and localized spins feel a quenched disorder. 

In the upper panel of Fig. [8] we show a comparison of the correlation function g{r) 
calculated in the presence of the diagonal disorder and without it. /^From this figure 
one sees that the influence of the diagonal disorder is almost negligible, at least for 
a small-to-moderate values of the potential V. For larger values of V, the presence 
of negatively charged out-of-plane oxygen atoms reduces the hole mobility resulting 
in a visible enhancement of the spin-spin correlation function. It is worthwhile to 
emphasize that the NDO constraint becomes very important in the presence of the 
diagonal disorder despite the low concentration of holes. In a homogeneous systems 
with (5 <C 1 this constraint is less important as the probability of double hole occupancy 
is proportional to S"^. Since the negatively charged out-of-plane oxygen atoms locally 
enhance the hole concentration, neglecting of the NDO may significantly modify the 
results [12]. 

Now we turn to the influence of the inhomogeneity-induced enhancement of the 
exchange interaction. Following the results of Ref. [12] we assume J to a be site- 
dependent quantity: 

J,, = J (1 + T],,) , (23) 

where 

{?7 > 0, if there is an out-of-plane oxygen atom 
above site i or j, 
0, otherwise. 

In contrast with the superconducting gap [42], the AF order is hardly modified by this 
mechanism. This can be clearly inferred from the lower panel in Fig. [HI The regime for 
a magnetic ordering predicted by many calculations in the t-J model extends to much 
larger dopings than observed in cuprates and this discrepancy is sometimes attributed 
to the inhomogeneities, which are neglected in many theoretical approaches (see the 
discussion in Ref. [E]). Although, inhomogeneities are expected to play an important 
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Figure 8. g(r) for J = 0.2t in tiie presence of the out-of-plane oxygen atoms. Tiie 
upper panel demonstrates the influence of the diagonal disorder, whereas the lower 
panel shows effects coming from the site-dependent exchange interaction. The model 
parameters (5, V, r/) are given in the legend. 



role in high-temperature superconductors [25], our resuhs indicate that their influence 
in the AF ordering is rather limited. In particular, we expect that the inhomogeneities 
introduced by the out-of-plane oxygen atoms cannot explain the above discrepancy. 

3.3. Finite-size effects 

Since our analysis has been carried out on finite clusters, it is necessary to check to 
what extent the results are affected by the finite size effects. One of the measures of the 
significance of the finite size effects is a sensitivity to the boundary conditions. Therefore, 
we have calculated the correlation functions and the spectral functions for systems with 
different boundary conditions and compared them to those which have been obtained 
with periodic boundary conditions. Here, we have used a method known as averaging 
over boundary conditions (ABC) [2H|- Each time a particular hole jumps out of the 
cluster, it is mapped back into the cluster with wave function with a different phase. 
Then the results are averaged over these phases; thereby the reciprocal space is probed 
in a much greater number of points than in the case of periodic boundary conditions. In 
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our simulations we have used a slightly modified version, where the phases were chosen 
randomly in each Monte Carlo step. The averaging over the boundary conditions has 
been carried out concurrently with averaging over the ensemble generated in Monte 
Carlo run [l3]. This way it does not require an additional computational effort. 

Another more direct way to check the infiuence of the finite size effects is to compare 
results obtained on clusters of different sizes. In order to do so, we have repeated 
some calculations on 40x40 cluster. Fig. [9] shows a comparison of spectral functions 
obtained for 20x20 and 40x40 clusters with periodic boundary conditions as well as 
for 20x20 cluster with ABC. Since the false-colour plots of these spectral functions are 
very similar to each other, we present their energy dependence for a few selected points 
of the Brillouin zone. One can see from this figure that the coherent part of the spectral 
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Figure 9. Spectral functions calculated for selected points of the Brillouin zone. These 
results have been obtained on 20x20 (solid line) and 40x40 (dashed line) clusters with 
periodic boundary conditions and on 20x20 cluster with ABC (dotted line). Since the 
positions of the coherent peaks obtained with help of these three approaches are almost 
indistinguishable, we have cut the vertical axis in such a way that the incoherent parts 
are more pronounced. These results have been obtained for J = 0.2t, kT = Q.lt, and 
S = 0.02. 



functions is almost exactly the same in these three cases. The low-intensity parts also 
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look very similar, thought some differences can be seen. Another quantity we have used 
to analyze the finite size effects is the correlation function g{r). Fig. [10] shows g{r) 
determined for the same three systems, for which the spectral functions are presented 
in Fig. [9l Despite minor quantitative differences the overall character of all three 

1.00 L I 




Figure 10. Correlation function g{r) for the same systems as in Fig. [Qj Also, 
the same parameters have been used. The insets show examples of snapshots of the 
spin configurations for 40 x 40 and 20 x 20 clusters [black (white) square corresponds to 

correlation functions is the same. Though we have not carried out a systematic finite 
size scaling, the similarity of both the spectral functions and the correlation functions 
constitutes a significant indication that our results are valid also in the thermodynamic 
limit. 

4. Summary 

We have developed a doped-carrier representation of the Ising t-J model. In this 
formulation, the system is described in terms of fermions interacting with static 
localized spins. Although it is a slave-particle approach, in contrast with many similar 
approaches, the local NDO constraint is taken into account exactly. The proposed 
Hamiltonian has the global Z2 symmetry at any values of the parameters, J and t. This 
model is of an interest in itself since it represents a simple though nontrivial electron 
system which captures the physics of strong electron correlations. The issue of how these 
correlations affect the magnetic ordering of the lattice spins is thouroughly investigated 
in the present work. Besides, this model may provide at least for some values of the 
model parameters a guess supported by unbiased numerical calculations regarding the 
actual low-energy behaviour of the quasiparticle excitations in a more realistic isotropic 
t — J model. 

In particular, we have calculated the one-hole energy and compared our results with 
those obtained for t-J^ and t-J models. We have found that the one-hole energy is the 
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same as in the t-J^ model in two limiting cases, J — and J — oo. For intermediate 
J the one-hole energy obtained for the t-J model is in between the results obtained 
for t-Jz and t — J\z2 approaches and the differences between t-J and t-Jz models are 
comparable to those between t-J and t — J\z2 ones. 

The main advantage of the present approach consists in that it can be applied for 
very large systems and the computational effort increases much slower with the size of 
the system than, e.g., for exact diagonalization. Moreover, it works for arbitrary value 
of the coupling J and for arbitrary doping level. In particular, in the small-J regime 
the exact diagonalization and Quantum Monte Carlo methods give rather poor results. 
This is because of the fact that in this regime the size of the defects generated by moving 
holes is comparable to or larger than the size of cluster the calculations can be carried 
on. Since the clusters in our approach are much larger, this problem is less significant. 
Additionally, our method does not require translational invariance of the system. This 
feature is especially important in the context of the recent experimental results, which 
clearly indicate the presence of inhomogeneities in curates. It could also be applicable 
to optical lattices, where the translational symmetry is broken by a trap. 

Using the proposed approach we have found that the AF SRO persists for 
temperatures and dopings which are much beyond the boundaries of the AF LRO phase, 
what is in agreement with recent experiments on the high-temperature superconductors. 
We have also demonstrated that the AF LRO depends on the exchange interaction J. 
It concerns the transition temperature as well as the maximal doping at which the AF 
LRO vanishes. We explicitly demonstrate that the local no double occupancy constraint 
plays the dominant role in destroying the magnetic order at finite doping. 

Finally, we have shown that the inhomogeneities induced by the out-of-plane 
oxygen atoms have a rather limited influence on the spin-spin correlations functions, 
at least in the underdoped regime. Although, localization of holes by their electrostatic 
potential stabilizes the AF LRO, this mechanism becomes important only for a relatively 
strong diagonal disorder. Such a limited influence of inhomogeneities on the AF order 
is closely related with the NDO constraint. Note that exactly at half filling the diagonal 
disorder does not influence the system, provided the NDO is properly taken into account. 
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